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Abstract 

In this paper, we study a non-linear numerical scheme arising from the im- 
plicit time discretization of the Bazant-Thonguthai model for hygro-thermal 
behavior of concrete at high temperatures. Existence and uniqueness of the 
time-discrete solution in two dimensions is established using the theory of 
pseudomonotone operators in Banach spaces. Next, the spatial discretiza- 
tion is accomplished by the conforming finite element method. An illustrative 
numerical example shows that the model reproduces well the rapid increase 
of pore pressure in wet concrete due to extreme heating. Such phenomenon 
is of particular interest for the safety assessment of concrete structures prone 
to thermally-induced spalling. 
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1. Introduction 



The hygro-thermal behavior of concrete plays a crucial role in the as- 
sessment of the reliability and lifetime of concrete structures. The heat 
and mass transfer processes become particularly important at high tem- 
peratures, where the increased pressure in pores may lead to catastrophic 
service failures. Since high-temperature experiments are very expensive, pre- 
dictive modeling of humidity migration and pore pressure developments can 
result in significant economic savings. The first mathematical models of con- 
crete exposed to temperatures exceeding 100°C were formulated by Bazant 
and Thonguthai in Since then, a considerable effort has been invested 
into detailed numerical simulations of concrete structures subject to high 
temperatures. However, much less attention has been given to the qualita- 
tive properties of the model, as well as of the related numerical methods. 

In particular, the only related work the authors are aware of is due to 
Dalik et al. |3], who analyzed the numerical solution of the Kiessl model [§[ 
for moisture and heat transfer in porous materials. They proved some exis- 
tence and regularity results and suggested an efficient numerical approach to 
the solution of the resulting system of highly non-liner equations. However, 
the Kiessl model is valid for limited temperature range only and as such it is 
inappropriate for high-temperature applications. In this contribution, we ex- 
tend the work [3j by proving the existence and uniqueness of an approximate 
solution for the Bazant-Thonguthai model, arising from the semi-implicit dis- 
cretization in time. A fully discrete algorithm is then obtained by standard 
finite element discretization and its performance is illustrated for a model 
problem of a concrete segment exposed to transient heating according to the 
standard ISO fire curve. Here, the focus is on the short-term pore pressure 
build up, which is decisive for the assessment of so-called thermal spalling 
during fire. 

At this point, it is fair to mention that the Bazant-Thonguthai model was 
later extended towards more detailed multi-phase description, see e.g. the 



works of Gawin et al. |6|], Tenchev et al. [13[ and Davie et al. |4| for specific 
examples. When compared to the original version, these advanced models 
provide better insight into physical and chemical processes in concrete (such 
as influence of gel water, pore water, capillary water, chemical reactions at 
elevated temperatures, etc.). Such potential increase in accuracy, however, 
comes at the expense of increased number of parameters, which typically 
reflect complex multi-scale nature of concrete. Hence, their experimental 
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determination is rather complicated and the parameters can often only be 
calibrated by sub-scale simulations. Therefore, in this work we adopt a prag- 
matic approach and consider the single-phase Bazant-Thonguthai model with 
parameters provided by heuristic relations, obtained from regression of reli- 
able macroscopic experiments. 

The paper is organized as follows. In Section 121 we present the gen- 
eral single-phase, purely macroscopic, model for prediction of hygro-thermal 
behavior of heated concrete. In Section [3j we introduce basic notation, the 
appropriate function spaces and formulate the problem in the strong and vari- 
ational sense. In Section H~T| we specify our assumptions on data and modify 
structure conditions to obtain a reasonably simple but still realistic model 
of hygro-thermal behavior of concrete at high temperatures due to Bazant 



and Thonguthai [1]. An application of the Rothe method of discretization in 
time leads to a coupled system of semilinear steady-state equations, which 
(together with the appropriate boundary conditions) form a semilinear ellip- 
tic boundary value problem, formulated in the form of operator equation in 
appropriate function spaces. The existence result for this problem in space 
W l,p (Q) 2 with p e (2,4) is proven in Section |4~21 using the general theory of 
coercive and pseudomonotone operators in Banach spaces. Next, the prob- 
lem is resolved using the finite element method as presented in Section 15.11 
In Section I5.2[ numerical experiments are performed to investigate the mois- 
ture migration, temperature distribution and pore pressure build up in the 
model of concrete specimen exposed to fire, including a simple engineering 
approach to study the spalling phenomenon. 

2. The coupled model for wet concrete 

2.1. Conservation laws 

The heat and mass transport in concrete is governed by the following 
system of conservation laws: 
energy conservation equation: 




-V ■ Jg(0,w,V6,Vw) + C w J. 



(8,w,V8,Vw)-V8; (1) 



dt 



w 



water content conservation equation: 
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The primary unknowns in the balance equations ([I])-© are the temperature 
9 and the water content w; w represents the mass of all evaporable water 
(free, i.e. not chemically bound) per m 3 of concrete. Further, H and M. 
represent the amount of (internal) energy and the amount of free water, 
respectively, in 1 m 3 of concrete, Jg is the heat flux, C w the isobaric heat 
capacity of bulk (liquid) water and J w the humidity flux. 

2.2. Constitutive relationships for heat and moisture flux 

Following [H, the heat flux Jg arises due to the temperature gradient 
(Fourier's law) and due to the water content gradient (Dufour flux) 

Jg(9, w, V0, Vw) = -D ew (9, w)Vw - D e0 (9, w)V9 (3) 

and the flux of humidity J w consists of the flux due to the humidity gradient 
(Fick's law) and due to the temperature gradient (Soret flux) 

J w {0, w, V0, Vw) = -D ww {9, w)Vw - D w9 {9, w)V9, (4) 

where D 9w , D 9e , D ww and D wd are continuous diffusion coefficient functions 
depending non- linearly on 9 and w. 

2.3. Boundary and initial conditions 

To complete the introduction of the model, let us specify the boundary 
and initial conditions on 9 and w. The humidity flux across the boundary is 
quantified by the Newton law: 

J w (9, w, V0, Vmj) ■ n = 7 c (w - w^), (5) 

where the right hand side represents the humidity dissipated into the sur- 
rounding medium with water content Woo, specified in terms of the film co- 
efficient 7 C . As for the heat flux, we shall distinguish the convective and 
radiation boundary conditions given by 

J e (9,w,V9,Vw)-n = ^(9-9^), (6) 
J e (9,w,V9,Vw)-n = a c (9 - 0«,) + ea(9\9\ 3 - 9 A J, (7) 

respectively, in which a c designates the film coefficient for the heat transfer, 
and 6*00 is temperature of the environment. The last expression in Eq. ([7j) 
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expresses the radiative contribution to the heat flux, quantified by the Stefan- 
Boltzmann law in terms of the relative surface emissivity e and the Stefan- 
Boltzmann constant a and the temperature difference (# 4 — 9^) Q The initial 
conditions are set as follows: 

0(0) = 9 , w(0) = w . (8) 

Here, 8q and wo represent the initial distributions of the primary unknowns 
6 and w, respectively. 

3. Notation and formulation of the problem 

Vectors, vector functions and operators acting on vector functions are de- 
noted by boldface letters. Throughout the paper, we will always use positive 
constants c, c\, C2, . . . , which are not specified and which may differ from 
line to line. For an arbitrary r G [1, +oo], L r (Q) denotes the usual Lebesgue 
space equipped with the norm || • ||L'-(n), and W k,p (fl), k > 0, p G [l,+oo], 
denotes the usual Sobolev space with the norm || • ||v^fc,p(n)- Let X be a Ba- 
nach space. By C([0,T],X) we denote the space of all continuous functions 
ip : [0, T] — > X. Throughout the paper p' = p/{p — 1), p > 1, denotes the 
conjugate exponent to p. <fi'{t) indicates the partial derivative with respect 
to time; we also write 4>'(t) = <p t . 

We consider a mixed initial-boundary value problem for a general model 
of the coupled heat and mass flow in a two-dimensional domain Q C M? with 
a Lipschitz boundary dQ, which consists of non-intersecting pieces and 
Tat, dfl = r^UT^r. r# re presents the part of the boundary which is exposed 
to fire, whereas the other part denoted by T N is exposed to atmosphere. Let 
T > be the fixed value of the time horizon, = O X (0,T), Tr? = 

x (0,T) and T^t = x (0, T). The strong formulation of our problem 



1 Replacing the term A with 6\6\ 3 is essential later in the proof of Theorem 2. 
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is as follows: 

H t = -V • J 9 + C W J W ■ V0 in Q T , (9) 

M t = -V ■ J w in Q T , (10) 

J w -n = j c (w-w 00 ) onrVrUTRT, (11) 

J e ■ n = a c (9 - Boo) on T NT , (12) 

J e n = a c (6 - + ea(|0| 3 - O on T RTl (13) 

0(0) = O in J], (14) 

w(0) = w in a (15) 

Here we assume that all functions are smooth enough. Now we can formulate 
the problem in the variational sense. Suppose that [0oo(t),iUoo(£)] £ C(0,T) 2 
and [0 O , w ] e W 1,r (fi) 2 , r > 2. Find a pair [0,w] G C([0, T]; W 1>r (fi) 2 ) such 
that 

- / (H, <f>'l(t)) + (M, dt - [ Jg ■ V0! + J w ■ V0 2 dQ T 

JO JQt 

- j C W J W ■ V00i dQ T + I ea(|0| 3 - 0^)0i dS T 

+ / Oi c {6 - tfoo)^! dS T + / lc{w - Woo)02 dS'r = 

JTrtUTnt JTrt^nt 

holds for all test functions (p = [&,0 2 ] e C °°([0, T]; C 00 ^) 2 ) and 

0(0) = O and w(0) = w in fi. (16) 

Such pair [9,w] is called the variational solution to the system (|9|)-( IT5|) . 

Remark 1. To the best of our knowledge, there are no existence results for 
the presented model available. 

4. Existence of the approximate solution to the Bazant's model 

4-1- Structural conditions and assumptions on physical parameters 
A\ In PJ, Bazant and Thonguthai expressed the time variation of amounts 
% and M, as follows: 

dU 36 dw d dw 

-m '= psCs m ~ hd ~dT - K W (17) 

dM dw dw d 

~bT '- "a ""aT (18) 
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Here p s and C s , respectively, are the mass density and the isobaric 
heat capacity of solid microstructure (excluding hydrate water), Wd 
represents the total mass of the free water released in the pores by 
drying. h a denotes the enthalpy of evaporation per unit mass and hd 
denotes the enthalpy of dehydration per unit mass. 

A 2 We assume that the parameters p s , C w , hd, a c , (3 C , a and e are real 
positive constants. 

A3 The functions C s = C s (9) and h a = h a (6) are positive continuous 
functions, Wd = Wd(0) is positive increasing function belonging to 
W /1,00 (R), 9 00 {t) and w^t) are given continuous functions of time and 
6 ,w eW 1 ' r (tt),r>2. 

A4 Water content w is connected with temperature T and pore pressure 
P via a so-called sorption isotherm w = $(^,P), which has to be 
determined experimentally for each type of concrete. We assume $ to 
be a continuous function such that $(£1, £2) — for £ e and $ = 
otherwise. 

A 5 Following jlj, we consider the cross effects to be negligible. This leads to 
simple phenomenological relations introduced by Bazant and Thonguthai 
in the form 

J e ■= -X C (9)V9 and J w := - ^^ VP, (19) 

where the thermal conductivity A c and permeability k are assumed 
to be positive continuous functions of their arguments and g is the 
gravitational acceleration. 

4-2. Solutions to the discretized problem 

Incorporating the relations (IT7j) - (ITUl) into the system (19"|)- (1T5]) we get 
the modified Bazant-Thonguthai model with primary unknowns w, 9 and P 
consisting of 
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conservation laws: 



dw „ fK(0,P)„„\ dw d (9) „ , x 

P.CW^ - M»)^ = V • (A c (9, P)Vff) 

_ C „^VF.V9 + fci ^> infe (21) 
5 dt 

state equation of pore water: 

w-®(P,9)=0 mQ T - (22) 
radiation boundary conditions: 

-A c (0, P)V0 • n = a c (0 - 0^) + ea(\6\ 3 6 - 6 4 J on T KT ; (23) 
Neumann boundary conditions: 

-A c (0, P) V0 • n = a c (0 - 0^) on rW, (24) 

_^ L P) V p. n = /3c(p _p oo) onr^UT^; (25) 

and initial conditions: 

P(0) = P and 0(0) = O in Q. (26) 

Let = t < ti < ■ ■ ■ < t N = T be an equidistant partitioning of 
time interval [0; T] with step At. Set a fixed integer n such that < n < 
N. In what follows we abbreviate <p(x,t n ) by ip n for any function if. The 
time discretization of the continuous model is accomplished through a semi- 
implicit difference scheme 

W n+ i - W n _ /«(fl 

lyP ^ + Wrf ( 6> "+ 1 ) ~ W d(°n) ^ 



At \ 9 + ) At 



p s C s (9 n ) 9n+1 At 9n - h a (9 n ) Wn+1 At Wn = V • (A c (0 n ,P„)V0 n+1 ) 



Cw <^ W P n • V0 n + h d Wd{d ^\- WM . (28) 
# At 
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Here, we assume that the functions 9 n , w n and P n are known. In what 
follows we study the problem of existence of the solution 9 n+ i, w n+ \ and P n +i- 
Incorporating the relation (122]) into the system f l20|) -f l26|) we can eliminate 
the unknown field w and consider the problem with only two unknowns 9 and 
P. Consequently, the existence of w n+ \ follows from the existence of 9 n+ i 
and -P n +i by the relation (|22|) . For the sake of simplicity we denote [71*, r] := 
[P n +i, 9 n+1 ]. Let us put k(x) = K,(9 n (x), P n (x))/g, X c (x) = X c (9 n (x), P n (x)), 
$ n = Q(P n ,9 n ) and introduce the functions 

Ri(x,k,t) = -^$(k,t) - ±w d (r), (29) 

1 11 

R 2 (x,tt,t) = —p s C s (9 n )T-—h d w d (r)-—h Q (9 n )^(iT,T), (30) 

Fi{x) = x® n ~~Kt Wd{9n)l (31) 

F 2 (x) = ^p s C s {9 n )9 n - ^w d {9 n ) - ^^U n - C w K{x)VP n ■ V9 n . 

(32) 

Obviously, we have to solve, successively for n = 0, . . . , iV — 1, the following 
semilinear system with primary unknowns [tt, t] 



-V • (k(x)Vtt) + R 1 (x,it,t) = F x {x) in Q, 
-V • (X c {x) Vt) + R 2 {x, 7T, r) = F 2 (x) in £1 



(33) 



and with the boundary conditions 

-A c (x)Vr • n = a c (r - 9^) + ea(\r\ 3 r - 9^ n ) on T R , 

-k{x)Vix ■ n = /3 c (tt - Poo in ) on T N U T R , 

-A c (x)Vr • n = a c (r - 9 oc „) on Y N . 



(34) 



Definition 1. The pair [k,t] G W 1,r (Q) 2 , r > 2, is called a variational 
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solution to the system (I3"3"]) - (1310) i_/f 

«:(a;)V7r- Vf 7r + A c (a;)Vr- Vf T da;+ / i?i(aj, 7T, r) iv + i^aj, 7r, r) v T dec 

+ / /3 c ir v n dS + / a c rv T dS + / ecr|r| 3 rf T dS' 
ir fl urjv ^r fl urjv </r fl 

= / F 1 (x)v 7T +F 2 (x)v T dx + / /3 c -Poo,n^7r + a c 6 l oo,nWrd5+ / ecr^ )n u T dS 
Jn Jr R ur N Jr R 

(35) 

holds for every [iv, i>tt] € VF 1 ' 7 " (fi) 2 , r' = r/(r — 1). 

The main result of this section is the following Theorem [2] and Corollary [2J 

Theorem 2. Assume that [P n ,9 n ] E W l ' p (VL) 2 with some fixed p E (2,4) is 
known and let A 2 -A 5 be satisfied. Then there exists the variational solution 
[7r,r] E W^iVt) 2 to the system (ES])- 



Proof. In order to prove Theorem [2], it is convenient to define the operator 
T : W^iO) 2 -»■ W~ l *{Sl) 2 by 

(T([tt,t]),[v w ,v t ]) = / («(aj)V7r) • Vv w dx + / Ri(x,tt,t) v n dx 

Jn Jn 

+ / (Xc(x)'Vt) ■ Vf T dx + / R 2 (x, it, t) u t das 
in Jn 



+ / p c Trv n dS + / a c rv T dS + / ecr|r| rv T dS 
Jr R ur N Jr R ur N Jv R 

(36) 

and the functional / E W~ l ' p (Vt) 2 by 

(/, [v^Vt]) = I F l (x)v n dx+ I F 2 (x)v T dx+ I ^P^v^dS 
Jn Jn Jr R ur N 

+ / (Xcd^nVrdS + / ea9^ n v T dS (37) 
Jr R ur N Jr R 

for all [v-n^Vr] E W 1,p '(ti) 2 . Since we assume [P n , 9 n ] E W l ' p (Q) 2 with some 
fixed p E (2,4), we have the following estimate for the convective term 

c w K(x)VP n ■ we n 

I (C w K,(x)VP n ■ W6 n )v T dx < Ci||P n || w -i lP ( n )||0 n || w -i lP ( n )||'y T || w i il ,/ (n) (38) 
Jn 
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for all v T G W 1,p '(ft). One can prove in the similar way that the other inte- 
grals in ( 13 7p are finite and the functional / is well-defined. The variational 
problem can now be treated as a single operator equation T([n,r]) = f. 
Obviously, / G W~ 1 ' p (fl) 2 implies / G W~ l ' 2 (Q) 2 . First of all we prove 
that for a given h G W~ 1,2 (ft) 2 there exists [n, r] G W l,2 (Vt) 2 : the solution of 
the equation C([tt,t]) = h with £ : W 1 ' 2 (fl) 2 W- 1 ' 2 (fl) 2 defined by (JM]) 
(substituting £ instead of T) for all [tv,tv] G W 1 ' 2 ^) 2 . 

Lemma 3. £ : W 1 ' 2 (fl) 2 -» W- 1 ' 2 (fl) 2 is bounded. 
Proof. Test §M3) by [tt,t] G W 1 ' 2 ^) 2 . Take into account A 2 -A A to get 
(£([7r,r]),[7r,r]) < c x |M|^i, 2(n) + c 2 ||r||^i,2 (n) 

L 2 (fl) + C 4|| r llL 2 (n) 

+/ 3 c||7r||i 2(a m + a c \\T\\ 2 L 2 m) + ea\\r ub 



(an) ' "cii ' 1 1 1,2 (an) tw||/ iiL 5 (dfi) 
< c 5|| [tt, ||^i,2 (q)2 +ecr||r||| 5 ( 9n ). 



Due to the trace theorem there exists a constant Q r such that 

IM|L9(sn) < ^11^11^1,2(0) for all v G W l ' 2 (ft), q>l. 
Hence C is bounded. □ 

Lemma 4. C : W 1,2 (f2) 2 — >■ W~ 1 ' 2 (fl) 2 is coercive. 
Proof. A 3 , A$ and the Young inequality yield 

#i(*,6,6)6 = (^$(6,6) -^(6)) 6 
= ^$(6,6)6-^(6)6 

> -r& ~ c(v) (^(6)) (39) 

for every 6,6 £ ^ an d arbitrary 77 > 0. Further, A 3 , A 4 and the Young 
inequality yield the existence of a positive function gi and a non-negative 
function g 2 (both of spatial variable x) such that 

i? 2 (*,6,6)6 = ^p s C s (^)^-^(^ d (6) + ^(^)$(6,6))6 

< gi (x)H -g 2 (x) \Jx g n, v[6,6] e K 2 . (40) 
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Now ([3SD, dSHD, (HOD, the embedding W 1 ' 2 ^) m> L 2 (fi) and the Friedrichs 
inequality imply 



(£([7r,r]),[7r,r]) = f k(x)\ Vvr| 2 da; + / /3 c |vr| 2 d5 



r R ur N 

X c (x)\Vr\ 2 dx + [ a c \r\ 2 dS + l ea\r\ 3 r 2 dS 
u Jr R ur N Jr R 

+ / R\(x, 7r, r) 7r da? + / i? 2 (x, 7r, r) r da? 
in in 

II ||2 I || ||2 || ||2 I || ||2 

^ c ilFllvK 1 ^(n) + c 2|| r llvy 1 . 2 (n) _ ^ll 7r llL 2 (n)+ c 3|l r llL2( f7 )-c 4 
> c 5 1| [7T,T] ||^i, a(n)2 -ce (41) 

with some positive constants ci, . . . , cq and choosing 77 sufficiently small. □ 
Lemma 5. C : W 1,2 {Vt) 2 — >■ iy~ 1,2 (f2) 2 zs pseudomonotone. 
Proof. Obviously, since k > and A c > in f2, the inequality 

^)(6-^) 2 + A c (a3)(6-e 2 ) 2 >0 
holds for all a; G and for all fa, &], e ^ ^ [&&]• □ 

Corollary 1. The smoothness assumptions on k, A c , i?i and i?2 (the smooth- 
ness of Ri and i?2 follows from A2-A4 and f l29|) and (130]) ) and Lemma [3H5] 
imply that £ is continuous, bounded, coercive and pseudomonotone. Now 
jlll Theorem 3.3.42] yields the existence of the solution [tt,t] E W l,2 (Vt) 2 to 
the equation C([tt,t]) = h for every h e W~ 1,2 (Q) 2 . 

To get higher regularity results go back and consider / e W ,p (tt) 2 C 
W /_1 ' 2 (^) 2 with some p G (2, 4) and rewrite the system (|33|) -( 134|) in the form 

-V ■ (7c(aj)V7r) = Fi(se) - R x {x,tt,t) in Q, 

-V^(A c (a;)Vr) = F 2 (x) - R 2 {x, vr, r) in fi, 

-A c (a;)Vr-n = a c (r - 6>oo )n ) + ea(|r| 3 r - J on T R , 
-k(cc) Vtt ■ n = /3 c (tt - P^ n ) on r K U T N , 

-A c (a?)Vr • n = a c (r - 6^) on IV 

(42) 
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It is easy to verify that for the functional / G W~ 1,p (ft) 2 C W _1 - 2 (fi) 2 
defined by (13"7j) and for the weak solution [ir, r] G VF 1,2 (0) 2 (whose existence 
is ensured by Corollary [1]) , the functional G G iy _1 ' p (f2) 2 given by 

(G, [v n , v T ]) = / (Fi - u w + (F 2 - P 2 ) w T da3 + / a c (r - fl^^dS 
./c Jr R ur N 

+ / P c (7T-P OQtn )v n dS+ [ (ea(\r\ 3 r-d^ n ))v T dS (43) 
./r fl urv Jr fl 

for every u T ] G W l,p ' (Vt) 2 is well defined. It is known (see j§, that 
for given G G W-^iyi) 2 defined by flU with p G (2,4) the Neumann 
problem for the elliptic system (|42p (where the right hand side is represented 
by G) possess the solution [n,r] G W 1,P (Q) 2 . This completes the proof of 
Theorem |2J 

Corollary 2. Following Theorem |2J [P & , fc ] G W 1 '^) 2 yields [P fc+ iA+i] G 
iy liP (fi) 2 with any p G (2,4). Since we suppose [P , O ] G W 1,r (fi) 2 with r > 
2, we can conclude, that [P n , n ] G M^ 1,p (f2) 2 successively for n = 0, . . . , iV — 1 
for any p G (2,r) if r < 4 and p G (2,4) if r > 4. Note that this solution 
needs not to be unique. 



5. Numerical results 

5.1. Finite element implementation 

Consider a polygonal approximation Q h to Q, defined by an admissible 
quadrilateral partition Q h = {Qi, Q%, . . . , Qn b } such that = U^ 1 Q e , 
Q h C Q and every element Q e has diameter at most 2h. N n is used to denote 
the number of nodes of the mesh and T^, or T 1 ^ stand for the part of the ap- 
proximate boundary T h = dVL h where the radiation and convection boundary 
conditions are prescribed. We associate with Q h a finite-dimensional space 
of piecewise bi-linear basis functions (recall that p G (2,4)) 

S h = G C° (^ k ^j '■ v \q £ ^ P2 and restriction to each edge of dQ e 

belongs to V x for e = 1, 2, . . . , iV e | C W l ' p (n h ), (44) 

where V s denotes the set of polynomials of degree < s, cf. 0, Section 5]. 
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From the implementation point of view, it is more convenient to derive the 
numerical scheme by considering all three unknowns (w,6,P) instead of the 
reduced version (135|) ; the approximate solution (w^ +1 , 9^ +1 , P^ +1 ) G [S h ] 3 is 
thus provided by the weak form of balance equations ( p7|) and (|28|) . tested 
by v w G S h and v P G S h , constrained by the isotherm relation A 4 enforced 
at the nodes. This leads to a system of non-linear algebraic equations 

■^C n (X n+ i — X n ) + K n X n+ i + R(X n+ i) = F n+ i, (45) 

where e.g. X„ +1 = (w n+1 , n+ i, P n +i) G ]R 3Ar ™ xl stores the unknown nodal 
values of water content, temperature and pore pressure at time t n +i, respec- 
tively. The constant matrices in (T4"5l) exhibit a block structure 

rwP \ / C w 

\ I r n+l 

K = | K ] ,F n+1 = ( F n+1 

(46) 

and the non-linear term reads as 

/ -R w (6 n+1 ) \ 
R(X n+1 ) = R e (0 n+1 ) . (47) 

\ w n+1 - $(0 n+ i, P„+i) / 

The sub-matrices in (146D and (T47j) are obtained by the assembly of element 
contributions, e.g. 

C=AC F: +1 =AF: +m , R*(0n + i)= A X(9 n+1 , e ). (48) 





























-, n,e' Ti+l -, n+l,e> 

e=l e=l 



e=l 



Here, A is the assembly operator 0, Section 2.8] and 9 n>e G R 3xl stores the 
temperature values at the nodes of the e-th element at time t n , related to 
the local approximation 



9* >e (x):=e»(x)\ Te = N e (x)0 nie , (49) 

in which N e : T e — > M lx3 denotes the operator of linear basis functions. 
Analogous relations hold for the remaining fields. 
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The individual symmetric positive-definite matrices C* G R 3x3 can now 
be expressed as 

= / N e (x) T N e (x)dx, (50a) 
C e Z = I h a (6 h n!e (x))N e (x) T N e (x)dx, (50b) 

J T e 

C = / (p s CMA X ))- h ^U X ))) N e( X ) T K(x)dx, (50C) 

whereas the symmetric positive-semidefinite blocks K* G IR 3x3 attain the 
form 

K P e = f ktW'^W)(VNeW) T VN e (,)dx 

JT e 9 

+ / (3 c N e (x) T N e (x)dS, (51a) 

Ke = l \\ C {eU X )i P U X ))^ e {x)) T VH e { X )dx 

JT e 9 

+ f a c N e (a;) T N e (a;)d5. (51b) 

J dT e n(r h N ur h R ) 

The non-linear terms R* G lR 3xl are expressed as 

= f w d(Ue(x)0)M e (x) T dx, (52a) 
R e e (9) = [ ea (N e (x)0) 4 N e (x) T dx (52b) 
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and the right hand-side-blocks F* e M 3xl are provided by 
K + x,e = ^ t J T W d (9 h nje (x))N e (x) r dx 



+ 



P c Poo(t n+l )N e (x) T dS, 



(53a) 



aT e n(r*ur*) 



r n+l,. 



a 



T e 9 



w k (6 h n /x),P* e (x)) Ve h n /x) ■ VP^ e (x)N e (x) T dx 



+ 
+ 



a^oo(t„+i)N e ( : r) T d5 



dT e n(r h N ur h a ) 

ea9Ut n+ i)M e (x) T dS. 



(53b) 



In the following example, the integrals were approximated using the 5x5- 
point Gauss quadrature and the non-linear system (jl5j) was solved iteratively 
using the Newton method with the residual tolerance set to 10 -8 . 



x *i 200.0 mm 
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Figure 1: Cross-section of specimen 



Figure 2: ISO fire curve 



5.2. Example 

Our model problem deals with a square concrete specimen of cross-section 
200 x 200 mm 2 , which is exposed to fire on the part Tr of the boundary and 
insulated on the remaining portion Tm (see Fig. [T}. We assume that the 
ambient temperature in the vicinity of Tr increases according to the ISO fire 
curve, 9oo(t) = 3451og(8t + 1) + 298.15, with t given in minutes (see Fig. |2}. 
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Physical quantity 


Notation 


Value 


Dimension 


Density of the solid microstructure 


Ps 


2400.0 


kg m~ 3 


Specific heat of liquid phase 


C w 


4181.0 


Jkg^K- 1 


Enthalpy of dehydration 


h d 


2.4 x 10 6 


J kg- 1 


Mass of cement per m 3 of concrete 


c 


300.0 


kg m~ 3 



Table 1: Material constants of concrete 



The uniform initial conditions are set to 9$ — 298.15 K and Pq = 2.7542 x 
10 3 Pa and the constants appearing in boundary conditions on r# are taken 
as a c = 25 WirtC -1 , (3 C = 0.019 ms" 1 , P M (t) = 2.7542 x 10 3 Pa, o = 
5.67 x 10" 8 Wm~ 2 K" 4 and e = 0.7. 

5.2.1. Material data for concrete at high temperatures 

Basic material constants for concrete employed in this example appear 
in Table [TJ Following j7j], the thermal conductivity of concrete A c can be 
bounded by lower A; and upper X u limit values defined by 

A,(0) = 2.0 - 0.2451((r3- 273.15)/100.0) + 0.0i07((f3- 273.15)/100.0) 2 , 
\ u {9) = 1.36 -0.136((^-273.15)/100.0) + 0.0057((^- 273.15)/100.0) 2 

and is set to \ c (9) = (Xi(9) + X u {9))/2 in the results reported below. 

1000 



100 / 

10 I 

J& RH=1.0 

1 7/ RH-0.8 

/ RH=0.6 

RH=04 

1 

298.15 398.15 498.15 598.15 698.15 
Temperature [K] 

Figure 3: Thermal conductivity of concrete Figure 4: Permeability of concrete 

The permeability of concrete k = k(9, P) is adopted from j7[, (12a)- 
(12b)] and the relative change in permeability k/kq [-}, k,q = 10~ 13 [ms -1 ], 
is displayed in Fig. H] as a function of temperature 9 and relative humidity 
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Figure 5: Thermal capacity of solid skeleton Figure 6: Mass of dehydrated water 



RH, denned as RH(P,9) = P/P sat (9), where P sa t{9) is the vapor saturation 
pressure (see ji])H 

The specific heat of solid matrix C s is considered in the form (cf. [H) 

C s {9) = 900.0 + 80.0(0 - 273.15)/120.0 - 4.0 ((6 - 273.15)/120.0) 2 . (54) 

The thermal capacity of solid skeleton is displayed in Fig. |5j Assuming 
that concrete is fully hydrated at room temperature, the mass of dehydrated 
water is given as [5( 

r 0.0 for 9 < 373.15 K, 

w d {9) = < 0Mc(9 - 373.15)/100.0 for 373.15 < 9 < 973.15 K, 
( 0.24c for 9 > 973.15 K, 

see also Fig. |6] for an illustration. The temperature dependence of the en- 
thalpy of evaporation h a will be approximated by the Watson equation 0| 



h a {9) = 2.672 x 10 5 (647.3 



\0.38 



(55) 



provided 9 < 647.3 K. Note that for higher temperatures there is no liquid 
water in the pores and h a {9) = 0. 

The last comment concerns the sorption isotherms. For relative humidi- 
ties RH < 0.96 and RH > 1.04, thermodynamics-based relations w = 
&(9, P) introduced in [jj are adopted. In the transition range, we employ 
the C 1 -continuous cubic interpolation. 



2 Note that the pore pressure P can generally exceed the saturation pressure P sa t, so 
that RH > 1. 
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Figure 7: Temperature, pore pressure and water content distribution 30 min (left) and 60 
min (right) after exposure to fire. 
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5.2.2. Results 

The results presented in this section are obtained using an in-house MAT- 
LAB code and correspond to the uniform spatial discretization by 80 x 80 
square elements (without any adaptivity) and to the time step At = 5 s. The 
distribution of individual fields at two characteristic times appears in Fig. [7J 

We observe that the numerical model correctly reproduces the rapid heat- 
ing of concrete specimen, accompanied by highly localized profile of water 
content distribution. Mainly the latter phenomenon, which is accurately re- 
solved by the adopted fine grid, then contributes to the development of high 
values of pore pressure in the interior of the structure, leading to a potential 
threat to its stability during fire due to explosive spalling. 




50 100 150 200 50 100 150 

x 1 [mm] x [mm] 



Figure 8: Spalling of concrete 40 min (left) and 60 min (right) after exposure to fire. Region 
(a) denotes elements failed due to spalling, (b) marks unstable elements and (c) stable part 
of the cross-section. 

In order to assess this behavior more quantitatively, we adopt a conser- 
vative engineering approach and assume that the spalling of concrete occurs 
when the (appropriately reduced) pore pressure exceeds the temperature- 
dependent tensile strength of concrete. In particular, the spalling at x e fl h 
and time t n occurs when j^, Section 3.3jf| 

4>P*{x) > f t (9 h n (x)) , (56) 



3 It should be emphasized that relation (|56|) is purely heuristic; it nevertheless corre- 
sponds surprisingly well with detailed numerical simulations cf. [12L Section 4.2]. 
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where is the porosity of concrete, P% and 9% denote the finite element 
approximations of pore pressure and temperature for mesh size h and the 
n-th time step and ft as a function of 9 is provided by a piecewise linear 
relation [5J, Section 4.4] 



where fto designates the initial tensile strength of concrete. 

The spatial distribution of the spalling damage for <p = 0.1 and fto = 
2 MPa appears in Fig. [8J Here, three different zones can be distinguished. 
The first region (a), highlighted in black color, corresponds to elements failed 
due to spalling damage as predicted by the criterion (|56|) for the element 
center. The second zone (b) corresponds to the part of the structure in 
which the local strength is sufficiently high to sustain the pore pressure, but 
its stability is lost due the explosive spalling of the former region. Finally, the 
light gray zone (c) indicates the portion of the cross-section still capable of 
transmitting stresses due to mechanical loading, which is thereby responsible 
for the structural safety during fire. 
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